NON- VARIATIONAL COMPUTATION 
OF THE EIGENSTATES OF DIRAC OPERATORS 
WITH RADIALLY SYMMETRIC POTENTIALS 



LYONELL BOULTON AND NABILE BOUSSAID 

Abstract. We discuss a novel strategy for computing the eigen- 
values and eigenfunctions of the relativistic Dirac operator with a 
radially symmetric potential. The virtues of this strategy lie on the 
fact that it avoids completely the phenomenon of spectral pollution 
and it always provides two-side estimates for the eigenvalues with 
explicit error bounds on both eigenvalues and eigenfunctions. We 
also discuss convergence rates of the method as well as illustrate 
our results with various numerical experiments. 



1. Introduction 

The free Dirac operator acting on 4-spinors of L^(]R^)'^ is determined 
by the first order differential expression 

3 

D := a ■ P + P = -i'^ atdt + /3, 
fc=i 

where a = [01,02,0(3), and the Pauh-Dirac matrices are: 

"^=(° 0) '^=(^0' -L 

forai=(° J)' ^^=(? 0') ^^=(0-1 

We assume that the units are fixed so that m = c = h = 1. Standard 
arguments involving the Fourier transform show that D defines a self- 
adjoint operator with domain if^(M^)^ and that the spectrum of D 
is 

Spec(L)) = (-CX), -1] U [1, 00). 

Spherically symmetric potentials are Hermitean 4x4 matrix multi- 
phcation operators, V, acting on /.^(M^)"^, such that C^(R^ \ {0}) C 
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Dom(V^) and 
where 



is the spin operator, and R is the matrix of the rotation of angle (p and 
axis n. Here Dom{V) denotes the maximal domain of V. 

Spherically symmetric potentials may be constructed from maps 
0sc,ci,am : ^ — > ^ via 

(1.1) V{x) = 0sc(k|)/? + 0ci(|a;|)/c4 + i0am(|a;|)/9a ■ ^. 

\x\ 

The subscripts "sc" , "el" and "am" , stand for "scalar" , "electric" and 
"magnetic" potential, respectively. Radial symmetry on the electric 
potential, for instance, is a consequence of the assumption that the 
atomic nucleus is pointwise and the electric forces are isotropic in an 
isotropic medium like the vacuum. In the particular case (psc = 0am = 
and (pciii^) = 1 /f-i \l\ < "\/3/2, H describes the motion of a relativistic 
electron in the field created by an atomic nucleus. 

If V is subject to suitable smallness and regularity conditions (such 
as those ensuring that V is relatively compact with respect to D), then 
H := D + V defines an essentially self-adjoint operator in C(f^(M'^\{0})^ 
with self-adjoint extension domain if^(]R^)^ and 

(1.2) Spec^jH) = Spec{D) = {-oo, -1] U [1, cx)), 

see |Tha92l Theorems 4.2 and 4.16]. In fact there are known conditions, 
satisfied by potentials of practical interest (see e.g. |BG87t Theorems 6 
and A]), preventing the existence of embedded eigenvalues. 

The addition of a non-zero potential might give rise to a non-empty 
discrete spectrum in the gap (—1,1). These eigenvalues can only be 
found explicitly for a few simple systems which do not go much be- 
yond the case of a single electron embedded in a field created by an 
atomic nucleus, see (14. ip . For more complicated potentials, one has to 
rely either on asymptotic techniques (cf. |BL94j . |GL99j . |Sch03j and 
references therein) or on numerical estimations. 

Few robust computational procedures are currently available to esti- 
mate numerically the eigenvalues of H, see |Dya90| , [DESOOj . |DESVOO] . 
|DES03] . |DEL07j . |LLT02] and references therein. Since H is strongly 
indefinite, a direct implementation of the projection method is not pos- 
sible due to variational collapse. Under no further restriction on the 
potential or the reduction basis, accumulation points of eigenvalues of 
the finite dimensional approximate operator which do not belong to 
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Spec(if) might appear, see |SH84] and |Dya90| . This phenomenon is 
known as spectral pollution. 

The aim of the present paper is to investigate the applicability of 
the so called quadratic projection method for finding eigenvalues and 
eigenfunctions of H. This method has been recently studied in an ab- 
stract setting (see [ShaOOj . |LS04j . |Bou06] and |Bou07j ) and it has 
already been applied with some success to crystalline Schrodinger op- 
erators, |BL07] , and magnetohydrodynamics, |Str08] . In this approach, 
explained at length in Section [31 the underlying discretised eigenvalue 
problem is quadratic in the spectral parameter (rather than linear), and 
has non-real eigenvalues. Its main advantage over a standard projec- 
tion method is its robustness: it never pollutes and it always provides 
a posteriori two-sided estimates of the error of computed eigenvalues. 

Section [3TT] is devoted to a self-contained description of the quadratic 
projection method. In Theorem [1] we present an alternative proof of 
Shargorodsky's non-pollution Theorem |LS04| Theorem 2.6]. Addition- 
ally we show that information about eigenfunctions can be recovered 
from the quadratic projection method, see 03.41) . 

In Section H] we test the practical applicability of the numerical 
scheme proposed in Section [3] by reporting on various numerical exper- 
iments. As benchmark potentials we have chosen the purely coulom- 
bic, sub-coulombic and inverse harmonic electric potentials. In order 
to perform these numerical experiments, we split the space into up- 
per and lower spinor component, after writing the problem in radial 
form. Spectral pollution in the standard projection method using this 
decomposition and the consequences of unbalancing the number of up- 
per/lower components is discussed in Section [2J 

We have chosen a basis of Hermite functions to reduce the continu- 
ous problem into finite-dimensional from. In Section 13.31 we perform a 
rigourous convergence analysis of the quadratic method using this basis. 
This analysis relies upon the general result |Bou07l Theorem 2.1]. The 
surprising numerical evidence included in Section 14.41 strongly suggest 
that, under appropriate circumstances, choosing a basis that heavily 
pollute the standard projection method, significantly improves conver- 
gence rates of the quadratic projection method. 

We have posted a fully functional Matlab code for assembling the 
matrices involved in the quadratic projection method for V a coulombic 
potential in the Manchester NLEVP collection of non-linear eigenvalue 
problems |BHM"'"08j . See the permanent link IBHM"*"] . In Appendix [Al 
we include details of the explicit calculation of all matrix coefficients. 
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2. Spectral pollution and upper/lower spinor component 

BALANCE 

Consider, for any \l/ G L^(R^)^, the spherical coordinates represen- 
tation: 

(2.1) ilj{r,6,(f)) = r\l/(r sin(^) sin(0), r sin(^) cos(0), r cos(6')) 

where {r,6,(f)) G (0, oo) x [0, vr) x [— 7r,7r). The map \1/ is an 

isomorphism between L^(]R^)^ and I/^((0, oo), dr) ® L'^{S'^)^. If we de- 
compose L'^{S'^Y the sum of the so called two-dimensional angular 
momentum subspaces ^m^.K^, the partial wave subspaces are given by 

^mj,K, = L^{{0, oo), dr) (g) Kmj,K,, 

so that L^(M^)^ = ^S^m ,K ■ Without further mention, here and below 
we always assume that the indices {mj,Hj) run over the set G 
{-J, ■ ■ ■ , j} and K, G {±(j + i)}, for j G : G N}. 

The r factor in (12.1 p renders a Dirichlet boundary condition at 0. 
The dense subspaces C(f (0, oo) ® C S)rnj,Kj are invariant under 

the action of H. If V is as in fll.ip . then if \ C^{0,oo) ® -^jAj 
unitary equivalent to 

.^^N ._ A+0sc(r)+0ei(r) -|: + ^ + 0am(r) \ 

^ "^•'^^ V f + V + <^am(r) -1 - 0sc(r) + Mr)) ' 

The operators Hmj,Kj are essentially self-adjoint in C^(0, oo)^ under 
suitable conditions on the potentials 0sc,ei,am- Then 

Spec (if) = |JSpec(ifm^,«.). 

Below we often suppress the sub-index {mj,Kj) from operators and 
spaces, and only write the index k = kj. Note that the eigenvalues of 
H are degenerate and their multiplicity is at least rrij. By virtue of 

Since Spec^gg{Hf,) = (— oo, —1] U [1, oo), are strongly indefinite. 

Let us consider a heuristic approach to the problem of spectral pol- 
lution for and the decomposition of L^(0, oo)^ into upper and lower 
spinor components. For simplicity we assume that the potential is 
purely electric and attractive, 0sc('^) = 0am('^) = and (f>ei{r) < 0. 

The pair (m, v) G Dom(iiK) is a wave function of if^ with associated 
eigenvalue E G (—1, 1) if and only if 
(2.3) 

{^^i+l-E)u+{-dr + -)v = and {dr + -)u+{(j)^i-l-E)v = 0. 
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System (12. 3p can be decoupled into 

(2.4) Leu = and v = -{(p^i - I - E)~^ (^dr + '^^ u, 

where 

Le := - [-dr + ^) (0el - 1 - [dr + ^) + (0el + 1 ' E) 



> (1 + E)-^ + + +0el + (l-^). 

If we assume (pci is relatively compact with respect to the non-negative 
operator {dr + ^)*{dr + ^), 

(2.5) min Spec,^^ Lij > {1 - E) > 0. 

Moreover, the expression for v in f l2.4l) yields u ^ and v ^ 0. Hence 
E G Spec^isc-^K if only if is in the discrete spectrum of L^;. 

Let A4n G L^(0, oo) be a nested family of finite-dimensional sub- 
spaces such that A4n C. J^n+i, 



[jMN = L\0,oo) 



N>1 

and Cnm ■= ® -Mri C Dom{Hf^). Let P/v denote the orthogonal 
projection onto M.n, so that — > / in the strong sense. With this 
notation we wish to apply the projection method to operator H,^ with 
test spaces Cnm- 

We first consider the case = M. Let {uiy,VNy G Cnn be a 
sequence normalised by||nAr|p + ||fAr|p = l, such that 

(2.6) PA,[(0el + 1 - En)un + i-dr + '^)vn] = 0, 

(2.7) P^lidr + ^)njv + (0ei - 1 - En)vn] = 

for a suitable sequence Ej^ E E (—1, 1). Let xn be such that 

Wat = -(0ei - 1 - -Eat)"-^ (^(9r. + ^) Mat + xa^. 
By virtue of CT . 

(2.8) P^(0el - 1 - ^^)XAr = 0. 

If we where able to prove that 

(2.9) ||PA.(-a, + -)xn\\ ^0, iV ^ cx), 

r 

by substituting into (12.61) . we would have ||P/v-^^£'MAr|| 0. Thus, by 
virtue of the min- max principle alongside with (12.51) . we would have 
G SpeCjisj,L^ and so G Spec^jg^ifK. 
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Unfortunately, (12.91) can not be easily verified. Since Pn — > / in 
the strong sense, (0ci — 1 — Ej^)xn — ^ in the weak sense. Since 
0oi < 0, xn also in the weak sense. Therefore we are certain 
that {—dr + j)xn and hence PnL^un — > weakly. This does 
not imply in general fl2.9p . it only gives indication that the latter is a 
sensible guess. 

The key idea behind the above heuristic argument motivates several 
pollution-free numerical methods for computing eigenvalues of if^, in- 
cluding the quadratic projection method discussed in the forthcoming 
section: as (12. 3p might be vulnerable to spectral pollution in the gap 
(—1,1) because (12. 9p is not necessarily guaranteed, we characterise the 
eigenvalues of this problem by testing whether is in the spectrum of 
an auxiliary operator. This auxiliary operator has essential spectrum in 
{Re(2;) > 0}, so that a neighbourhood of is always protected against 
spectral pollution (see the discussion preceding [Bou06[ Lemma 5].) In 
the quadratic projection method, (12.91) is substituted by (13.40 . 

Remark 1. A successful procedure for computing eigenvalues of the 
Dirac operator is the one developed by Dolbeault, Esteban and Sere, 
[DESOOj and |DES03j . This procedure systematically implements the 
above idea. Multiplying by u the left side equation of (12. 4p and inte- 
grating in the space variable, gives A{E)u = for 

A{X)v:= +(0ei + l-A)K;|Mr. 

Both terms inside the integral decrease in A so, if v is regular enough, 
there is a unique A = \{y) G M satisfying A{\)v = 0. Upper estimates 
for the eigenvalues of if^ are then found from those of the matrix 
corresponding to the A-dependant form A{X)v restricted to f G A4n- 

Suppose now that M : N — > N. Let (n7v,fAr)* G C^MiN) be a 
sequence normahsed by||nAr|p + ||i;Ar|p = l, such that (12. 6p holds true 
and Pn is replaced in ([221) by Pm{n)- If limN^^ M{N)/N < 1, we 
are certainly less likely to obtain (12. 9p as P/v is replaced in (12.81) by 
Pm(N)- If, on the other hand, \mi]si-^c>o M[N) / N > 1, then we would 
be more confident about obtaining (12.90 for symmetric reasons. In 
Section we will present a series of numerical experiments supporting 
this heuristic argumentation. In particular, spectral pollution in the 
standard projection method appears to increase as ImiN^ao M{N) /N 
decreases (see figures H]) . 
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3. The quadratic projection method 

3.1. The second order relative spectrum. We now describe tlie 
basics of the quadratic projection method. In order to simplify the 
notation, below and elsewhere G denotes a generic self-adjoint operator 
with domain Dom(G) acting on a Hilbert space. One should think of 
G as being any of the if^ introduced in the previous section. The inner 
product in this Hilbert space is denoted by (■, ■) and the norm by || ■ ||. 

Let £ C Dom(G) be a subspace of finite dimension. Assume that 
C = Spanjfei, . . . ,bn}, where the vectors bj are linearly independent. 
Let 

K:= {{Gb,,Gb,))l,^,, L := {{Gb, , b,))lk=i 
and 5:=((6„M)-fc=i- 

For 2 e C, let Q{z) := Bz"^ - 2zL + K G C"^". The aim of the 
quadratic projection method is to compute the so called second order 
spectrum of G relative to C: 

Spec2(G, C) := Spec(Q) = {A G C : Q{X)v = 0, some O^ve C"}. 

Since i? is a non-singular matrix and det Q{z) is a polynomial of degree 
2n, Spec2{G,C) consists of at most 2n points. These points do not lie 
on the real line, except if C contains eigenvectors of G. However, since 
Q{zy=Q(z), 

Spec2(G',£) = Spec2(G, £). 

Approximation of the discrete spectrum of G using the second order 
spectrum has been discussed in |LS04j . |Bou07] . |BL07j and the refer- 
ences therein. 

The following result establishes a crucial connection between Spec(G') 
and Spec2(G, i^). Without further mention, we often identify the ele- 
ments V E C with corresponding v^C^in the obvious manner: 

n 

v = Y,{^,b*)b, and v = {{v,b*))]^„ 

k=l 

where {b*} is the basis conjugate to {bj}. Note that if the bj are 
mutually orthogonal and \\bj\\ = 1, then b* = bj. Below and elsewhere 
Us denotes the orthogonal projection onto a subspace S C Dom(G'). 

Theorem 1. Let C C Dom(G') be any finite- dimensional subspace. If 
A G Spec2{G,C), then 

(3.2) [Re(A) - |lm(A)|,Re(A) + | lm(A)|] n Spec(G) ^ 0. 
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Moreover, suppose that E is an isolated eigenvalue ofG with associated 
eigenspace S C Dom(G'). Let 

dE '■= dist(ii^, Spec G \ {E}) = min{|£' — x\ : x G Spec G, x E}. 

If 

(3.3) [Re(A) - |Im(A)|,Re(A) + | Im(A)|] n Spec(G) = {E} 

and Q{X)v = for ^ vE C", then the corresponding v E C satisfies 

(3.4) \\v-Ilevl^\lm\\ 



.\v\\ dE 

Proof. Let A G Spec2{G,C) and assume that Im(A) ^ 0. Since 

Q{z) = {{{z-G)b„{z-G)b,))l,^„ 
Q{X)v_ = for non-trivial w G if and only if 

(3.5) {{X-G)v,(X-G)w) = ^weC. 
For u,w & C and z = fi + iu where fi,!/ ^M., we have 
{{z—G)u,{z—G)w) = {{ji—G)u, {jji~G)w)—u^{u,w)—2iu{{jji—G)u,w). 
In particular, if we take w = v m (13.51) . we achieve 

||(ReA - G)t;f - | Im AnK;f - 2i\ ImA|((Re A - G)v,v) = 0. 
Thus 

(3.6) M5ilz^ = |„„,| 

and 

(3.7) |ImA|((ReA-G')?;,t;) = 0. 
But recall that G = G*, so 

{x-G)u\\ 



dist(x, Spec G) = min 



«GDom{G) 

for X G M. Therefore 

dist(Re A, Spec G) < |IniA|, 

confirming (13.21) . 

For the second part, assume that A, E and S are as in the hypothesis, 
and let v satisfy (13.51) and hence (13.61) . Then 

\\{E-G)v\\ < |E- Re Allien + ||(ReA-G)v|| 
< 2|ImA|||i;||. 
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Since 

\\{x-G)u\ 



we have 



so that 



dist(x, SpecG \ {E}) = min 

ueDom{H),u±S \\U\ 



\\{E-G){I-Ue)v\\>dE\\{I-Ile)v\\, 
\\iE-G)v\\ 



ll(/-n,)^||< , 

dE 

Now, 

\\{E-G)vf = \E-ReX\^\\vf+\\{ReX-G)vf+2{E-ReX){{ReX-G)v,v). 
Thus, by ([32]) and (lO) . 

||(^ - G)vf = \E- ReX\^\\vf + I lmX\^\\vf = \E - X\^\\vf. 



This gives 



11 / ^ ^ N ,, I A| ,, , 



dL 

as needed. □ 

This theorem suggest a method for estimating Spec(G) from 
Spec2{G,C). We call this method the quadratic projection method. 
Choose a suitable C C Dom(G), find Q{z) and compute Spec2(G',£). 
Those A G Spec2{G, C) which are close to M will necessarily be close 
to Spec(G), with a two-sided error given by | Im(A)|. Moreover if A is 
close enough to an isolated eigenvalue E of G, then 

(3.8) |ReA-E| < |lmA| 

and a vector ^ v E C such that Q{X)v = approaches the eigenspace 
associated to this eigenvalue with an error also determined by | Im(A)|. 
Note that there is no concern with the position of E relative to the 
essential spectrum, or any semi-definitness condition imposed on G. 
The procedure is always free from spectral pollution. 

Remark 2. A stronger statement implying the first part of Theorem [1] 
can be found in jLS04l Lemma 5.2]. In fact, for isolated points of the 
spectrum, the residual on the right of (13. 8p can actually be improved to 
for I Im(A)| sufficiently small, see |BL07| Corollary 2.6]. How- 
ever, note that this later estimate is less robust in the sense that is 
not known a priori. 
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3.2. The Hermite basis. In the forthcoming sections we apply the 
quadratic projection method to G = iJ^- To this end we construct 
finite-dimensional subspaces L C Dom(i?^) generated by Hermite func- 
tions. 

Let the odd-order Hermite functions be defined by 

-1 -ri 
^k{r) := C2fc+i/i2fc+i(r)e 2 , r > 0, 

where hnir) are the Hermite polynomials and c„ = \j2'^~^n\^ are nor- 
malisation constants. Motivated by the results of [Bou07l Section 3.4] 
on Schrodinger operators with band gap essential spectrum, we choose 

Below we might consider an unbalance between the number of basis 
elements in the first and second component, 7^ M. Without further 
mention, we often write Cn = (^nn, and £„, = CN{n),M{n) when M and 
depend upon n. 

We now compute the matrix coefficients of Q{z). For this we recall 
some properties of $A;(r). The Hermite polynomials are defined by the 
identity 

Jn 

K{z):={-ire^"—e-'\ n G N. 
They satisfy the recursive formulae, 

(3.10) h'Sz)=2nK_^{z), 

(3.11) K+i{z) = 2zhn{z) - 2nhn-i{z) 

and form an orthogonal set on the interval (—00, 00) with weight factor 
e , 

/ hm{z)hn{z)e ^ dz = T'n\^bnm- 

J —00 

The generating function of this family of polynomials is 

00 

^K{z)'-_ 



j-n 



n=0 



Thus 
(3.12) 



^ r ^ (-D'^i^n n-odd, 

>/ln(0)-r=>^ , SO that /l„(0)=^ 



n=0 n=0 
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The odd-order Hermite functions are the normahsed wave functions 
of a harmonic oscillator, 

(3.13) - $'^(r) + r^^k{r) = {2k + l)$fc(r), 

subject to Dirichlet boundary condition at the origin. They form an 
orthonormal basis of L^(0, oo) and so = J in (13.11) . 

The entries of the matrices K and L in (13. ip can also be found ex- 
plicitly from known properties of the Hermite polynomials. The crucial 
terms for assembling these matrices are obtained in Appendix \^ 

3.3. Convergence. The procedure described above is useful, provided 
we can find points of Spec2 (-&«;, C,) near the real axis. Here we formu- 
late sufficient conditions on the sequence of subspaces in order to 
guarantee the existence of a sequence A„ G Spec2(-ffK5 -^n) accumulat- 
ing at points of the discrete spectrum of -ff^. We then show that the 
sequence of subspaces (13.91) satisfies these conditions. 

We firstly recall the following result [BouOTj Theorem 2.1]. Note 
that Cnm are subspaces of Dom(if^) for all k. 

Theorem 2. Let E he an isolated eigenvalue of finite multiplicity of G 
with associated eigenspace denoted by S. Suppose that Cn C Dom(G^) 
is a sequence of subspaces such that 

(3.14) \\GP{u-Unu)\\ <6{n)\\u\\ Vu G ^, p = 0, 1, 2, 

where 6{n) ^ as n ^ oo is independent ofu and p. Then there exists 
b> and Xn G Spec2{G,Cn), such that 

(3.15) \Xn - E\ < b6{ny/\ 

We now verify (13.141) for G = if^ and as in (13.91) . To this end 
we consider an argument similar to the one discussed in [Bou07t Sec- 
tion 3.4] for the case of the cristaline non-relativistic Schrodinger op- 
erator. 

Let 

'A 0' 




A 



where A = —d^ +r'^ acting on L^(0, oo), subject to Dirichlet boundary 
conditions at the origin. We fix the domain of A as 

Dom(A) = U Cn, 

JVgN 

so that A = A*, A has a compact resolvent and Cn Q Cn-i are the 
eigenspaces of A. 



12 LYONELL BOULTON AND NABILE BOUSSAID 

For / and g regular enough, we have 

Hjf] = f^i/ + M and Hi f ^"l = + 
\9 J \Pzf + Pi9) ''\9j yq-if + Qi9 

where pj are linear polynomials and qj are quadratic polynomials in the 
variables dr, 0sc,ei,am and k/v. Condition fl3.14p is achieved by showing 
that both and are relatively bounded in the sense of operators 
with respect to A. 

The following results can be easily extended to more general poten- 
tials. Here we only consider those of interest in our present discussion. 

Lemma 3. Suppose that 
(3.16) 

J, _ , / h S IXsc,ci,am(r) I < cr \ Vr > 0, 

Vsc,cl,am Xsc,cl,am ~r 'V^sc,el,am WiieTe ^ „;, ^ roo/n 

[ "V^sc.cl.am t (,U, OOj, 

for some constant c > 0. Then Dom(y4) C Dom(if^) and there exist 
constants a,b > such that 

\\HPv\\ < a\\Av\\ + b\\v\\ E Dom(v4), p = 0,l, 2. 

Proof. It is enough to check that the Hj! are relatively bounded with 
respect to —d^Ic^. This, on the other hand, is a straightforward con- 
sequence of Hardy's inequality. □ 

By combining this lemma with Theorem[2]and Theorem[T], we achieve 
the following. 

Corollary 4. Let E he an isolated eigenvalue of H^^ of finite multiplic- 
ity with associated eigenspace E . Let L^u be defined by 03. 9p . Suppose 
that 0sc,ei,am Satisfy (13.161) and assume additionally that E C Dom(A^) 
for some q > 1. For b > large enough and independent of N or M, 
we can always find a sequence \nm £ '^'9QC2{H^, Cnm) , such that 

\Xnm - E\ <b{N~'^ + M~^) and 



(3.17) 

where vnm G Cmn solves Q{Xnm)vnm = 0. 



II (/ - H^)t;7VA/|| < biN~"-^ + M-"^; 



Proof. Let v = ^ ^ normalised by ||t;|| = 1. By virtue of 

Lemma [3] and Theorem [2], the desired conclusion follows if we are able 
to find 6 > such that 

oo 

(3.18) J]|(i^'i;„$,)|2<6n-2{«-'), 

k=n 
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for r = 0, 1 and j = 1,2. In order to show f l3.18p . note that the 
hypothesis v G Dom(y4'^) ensures 

oo oo 

k=n k=n 
oo 

k=n 

as n — > cxD. □ 

Remark 3. Suppose that the number of degrees of freedom M + N = 
2n is fixed. Modulo the constant b, the bound on the right side of (13.171) 
is optimal at {N,M) = {n,n). This suggests that an optimal rate 
of approximation might be achieved by choosing an equal number of 
upper/lower spinor components in fl3.9p . Contrary to this presumption, 
and depending on the potential V, the numerical evidence we present 
in sections 14.41 and [4.51 show that the residual on the right side of (13.81) 
can in some cases decrease significantly (over 18% in some cases) by 
suitably choosing 7^ M. 

We now explore precise conditions on the potential, in order to guar- 
antee the hypothesis of Corollary HI 

Lemma 5. Let (/'scci.am ^ C°°(0, 00) be such that 0sc,ei,am('^) ^ as 
r 00. Assume additionally thatr 1— > r'^4'sc,c\,am{i^) ore locally bounded 
for some a G (0, 1). Let H^u = Eu. For sufficiently small a > 0, 

||e'"^n||HP(o,oo) < 00 Vp G N. 

Proof. See |BG87l Corollary 3.1]. Let 

For any < 5 < min | ± 1 1 , we can always separate V = + 
where: 

(a) Vc is smooth, it has compact support and a singularity of order 
0(r^") at the origin, 

(b) Ve is smooth with bounded derivatives and ||V^||oo < 
Then 

u = -{D^ + Ve-E)-^V,u. 
Multiplying this identity by -D^ and e"^, yields 
(3.19) Dle^^^u = -DP{D^ + V, + ta- E)-^e'''^V,u. 
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Let o" G M. On the one hand, (D^ + Vs + ia — E)^^ is a bounded 
operator from H"{0, oo) to if'^"'~^(0, oo). On the other hand, multiph- 
cation by Vc is a bounded operator from i7°"(0, oo) into i7°"~"(0, oo). 
Indeed, note that [Tri99l Theorem 2.1(i) for p = q = 2] yield the latter 
result for cr G N after commutation and iteration, then duality and in- 
terpolation ensure it for all a G M. Therefore, as {DK + Vs + ia — E)~^Vc 
is bounded from iJ'"(0, oo) to iJ'^+(i~°)(0, oo), (KT^ and a standard 
bootstrap argument ensure the desired conclusion. □ 

We remark that, by using [HisOOl Lemma 5.1] (which can be modified 
in order to get ride of the term | in its assumption (i)), necessarily 

a < a/(1 — sy — E'^ for the above lemma to hold true. 

Corollary 6. Assume that 0sc,ei,am o,re as in Lemma\^ Let E he an iso- 
lated eigenvalue of Hf, of finite multiplicity with associated eigenspace 
S. Let Cnm be defined by (13. 9p . For b > large enough and indepen- 
dent of N or M, we can always find a sequence \nm ^ Spec2(-ffK, jCnm), 
such that (I3.17P holds true if q < 5/4. 

Proof. Let f G £^ be as in the proof of Corollary |H We show that 

oo 

(3.20) J2('^k + l)'^{v„<^,)\' <oo, 

k=0 

for j = 1, 2. Let F = {dr — r). Integration by parts and (I3.10p . yield 
1 

{vj,(^k) = / Vj{r)h2k+iir)e-'''2dr 

C2fc+1 Jo 

1 , _rl 

Vj{r)h'-2k+2{r)e ^ dr 



2{2k + 2)c2k+i 
-1 

2(2A; + 2)c2fc+i 
1 



22(2fc + 2)(2A; + 3)c2fc+i 
1 

22(2A; + 2)(2A; + 3)c2fc+i 
-1 



/ Fvj{r)h2k+2{r)e 2 dr 
Jo 

/ Fvj{r)h'2k+3{r)e 2 dr 
Jo 

i F Vj{r)h2k+?,{r)e 2 dr 
Jo 

POO 2 

F^v,{Q)h2k+A{^) + / F^Vj{r)h2k+A{r)e-'^ dr 
Jo 



23(2A; + 2)(2A; + 3)(2A; + 4)c2A:+i 



F^v,{Q)h2k-,m ^ /o°°FS(r)W(r)e-Vdr 

+ -T J—-. TT — «! + a2. 



2\2k + 2)(2A; + 3)(2A; + 4)c2fc+i 24(2A; + 2) . . . (2A; + 5)c2fc+i 

Identity (13.120 alongside with the fact that \F'^Vj{Q)\ < 00 (see Lemma[5]) 
and the Stirling formula, ensures that ai ~ n"''/^ as n ^ 00. On the 



NON- VARIATIONAL COMPUTATION OF EIGENSTATES 



15 



other hand, Lemma [5] ensures that F^Vj G L^(0, oo). Thus, since 

C2fc+5 -2 

~ n 

2\2k + 2)...{2k + 5)c2fc+i 

a2 = 0(?7,~^) as n ^ oo. This guarantees fl3.2UI) for g < 5/4. □ 

According to this corollary, for smooth potentials, the order of ap- 
proximation of the quadratic projection method to any eigenvalue E 
of should be at least a power 1/8 of the dimension of Cmn- The 
numerical experiment performed in Section 13.31 show that this bound 
improves substantially for particular potentials. See Figure [H] (right) 
and Table [5] (right). 

Remark 4. The arguments involved in the proof of this theorem show 
that, the behaviour of the wave functions at the singularity and its 
regularity are the main ingredients responsible for controlling the speed 
of approximation when using a Hermite basis fl3.9p . 

4. Some numerical experiments 

We now report on various numerical experiments performed for very 
simple radially symmetric potentials. It is not our intention to show 
accurate computations, but rather to illustrate how the method dis- 
cussed in Section 12] can be implemented in order to rigourously enclose 
eigenvalues and compute eigenf unctions of the Dirac operator. 

4.1. Ground state of the purely coulombic potential. We begin 
by considering the analytically solvable case of V being a radially sym- 
metric purely coulombic potential: 0sc = 0am = 0, (pdij) = 7/r. Here 
— V3/2 < 7 < 0. The eigenvalues of H,^ are given explicitly by 

/ 2 \ 

(4.1) E,= il + 7 

Note that Ej ^ 1 as j ^ cxd for all values of k. The ground state of the 
full coulombic Dirac operator H is achieved when k = — 1 and j = 0. 

In Figured] we superimpose the computation of Spec2{H^i, Cn) for 
two values of n, in a narrow box near the interval [—3,3]. Here £„ 
is given by (13.91) . For the set of parameters considered (k = — 1 
and 7 = -1/2), (jH]) yields Eq ^ 0.866025, Ei ^ 0.965925, E2 ^ 
0.9851210, E3 ^ 0.99174012 and E^ ^ 0.9947623. A two sided ap- 
proximation of Eq is achieved from the point A G Spec2(-ff-i, >Ciooo) at 
A ~ 0.8661 -|-0.0236i. According to (13.21) . there should be an eigenvalue 
ofH_i in the interval [0.8661-0.0236,0.8661 + 0.0236]. This eigenvalue 
happens to be Eq. For Ei, E2 and the pair (£'3, E4), we can also derive 
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similar conclusions. Note that SpeC(,gg(iir_i) is also revealed by points 
of Spec2(-f^_i, C-n) seemingly accumulating at (—00, —1] U [1, 00). 

In Figure [2] we show approximation of the corresponding ground 
wave function associated with Eq. We have also depicted the analytical 
eigenfunction: 



where is chosen so that ||mo|| = 1. From this picture it is clear that, 
at least qualitatively, UQ{r) seems to be captured quite well even for 
small values of n. 

We show a quantitative analysis of the calculation of Uq in Table O 
In the middle column we compute the residual on the left of (13.41) and 
on the last column we compute the right hand side of (13. 4p . It is quite 
remarkable that the actual residuals are over 74% smaller than the 
error predicted by Theorem [H 

4.2. /3-Dependence of the sub-coulombic potential. We now in- 
vestigate the case of the potential being radially symmetric and sub- 
coulombic: (psc = 0am = 0, (pci{r) = for P e (0,1). Here 
— 1 < 7 < 0. The purpose of this experiment is to show how The- 
orem [1] provides a priori information about Spec^^^^^Hf.) even for small 
values of n. Note that (11.21) is guaranteed from |Tha92t Theorems 4.7 
and 4.17]. Furthermore H has infinitely many eigenvalues according to 
|Tha92l Theorem 4.23]. 

In Figure [31 we show computation of the ground state of H_i, for 
p = 0.1: 0.1 : 1 and 7 = -1/2. As /3 ^ 1", Eq 0.8661, the ground 
eigenvalue of the coulombic Dirac operator. As (3 the eigenvalue 

remains above 1/2. Note that the family of operators is not analytic 
at /3 = for this potential. For (3 = the spectrum becomes 



The vertical bars show | Im(A)|, the maximum error in the computation 
of Eq ^ Re(A) given by Theorem [TJ For this example we have chosen 
n = 15. Table H] contains the data depicted in Figure [2l Observe 
that the error increases as /3 ^ O"*" and /3 — 1~. This seems to be a 
consequence of the fact that Eq becomes closer to other spectral points 
at both limits, so (Ieq ^ 0. 

4.3. The inverse harmonic electric potential. In this set of exper- 
iments we consider another canonical example in the theory of Dirac 
operators: 0sc = (f>a.m = 0, 4>ei{r) = + r^) for 7 < 0. The discrete 
spectrum of H is known to be finite for —1/8 < 7 < and infinite 



(4.2) 




Spec(i^«) = (-00, -3/2] U [1/2, 00). 
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for 7 < —1/8, [KlaSOj . As the parameter 7 decreases, we expect that 
eigenvalues will appear at the threshold 1, move through the gap, and 
leave it at -1. This dynamics is shown in Figure HI for the ground eigen- 
value of H_i. It is a long standing question whether the eigenvalues 
become resonances when they re-enter the spectrum. 

In Figure [5] we depict the first three eigenfunctions of H^i for 7 = 
—4. They correspond to the eigenvalues Eq ^ —0.3955, Ei ^ 0.6049 
and E2 ^ 0.9328. See also Figure H Note that for 7 = -4, both 
components of the eigenfunctions appear to obey a Sturm-Liouville 
type oscillation hierarchy. 

4.4. Upper/lower spinor component balance and approxima- 
tion of eigenvalues. We now investigate the effects of "unbalancing" 
the basis by choosing N ^ M. 

In Figure O we have performed the following experiment. Fix the 
number of degrees of freedom, dim{CMN) = 200. Then for = 10 : 
5 : 190 and M = 200 — N, use the quadratic method as well as the 
Galerkin method to approximate eigenvalues of in the spectral gap 
(— 1, 1). We firstly consider (^sc = 4>a.m = and (l)ei{r) = — l/(2r). 

The Galerkin method might or might not produce spurious eigenval- 
ues. The quadratic method will always provide two-sided non-polluted 
bounds for the true eigenvalues with a residual, obtained from (13. 2p . 
which might change with A^. See also figures |4] and [61 The Galerkin 
method appears to pollute heavily near the upper end of the gap for 
A^ > M, as predicted by the considerations of Section [2l Moreover, 
for the ground state, the minimal | lm(A)| is not achieved at A^ = 100 
which corresponds to A^ = M, but rather at some A^ > 100. It is 
remarkable that the residual are reduced significantly (up to 66% for 
the true residual) when M{N)/N ^1/5. 

If we performed the analogous experiment for the inverse harmonic 
potential, the conclusion are also rather surprising. See Figure [71 The 
Galerkin method appears to pollute heavily near the upper end of the 
gap for A^ > M as predicted in Section [2l However, now the ap- 
proximation is improved by over 16% for Eq and over 18% for Ei, if 
M{N)/N ^ 3. 

We can explain these phenomena by considering the relation between 
the components of the exact eigenvectors. 

In the case of a purely coulombic potential, the ground state is given 
by (14.21) where uo is a real constant. The lower spinor component just 
differs from the upper one by a scalar factor. When 7 G (0,1), the 
lower component is smaller in modulus than the upper one. Choosing 
A^ > M, can reduce an upper bound of the residual associated to the 
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first component, while the residual associated to the second component 
remains small due to the smallness of the lower component. 

In the case of an inverse harmonic purely electric potential, this argu- 
ment fails, as the two spinor components of the eigenfunction are not a 
scalar factor of each other, see Figure [51 If we denote an eigenfunction 
hju = (m"pp,m'°"), the figure suggests that |<92m"pp(0)| « |(92m1°™(0)|. 
As |(9^$fc(0)| = 0, it is natural to expect that a decrease in the residual 
is only achieved by choosing a suitable M > N. 

Remark 5. Although we can not prove it rigourously, strong evidence 
suggests that (for any of the potentials considered above) no spuri- 
ous eigenvalue is produced by the Galerkin method when N = M. 
Why bothering then with more complicated procedures, such as the 
quadratic projection method, to avoid inexistent spectral pollution. A 
partial answer is, on the one hand, robustness: we do not know a priori 
whether the Galerkin method pollutes for a given basis. On the other 
hand, as the experiments of this section suggest, some times forcing a 
kinetic unbalance into a model might improve convergence properties. 

4.5. Convergence properties of the odd Hermite basis. A con- 
vergence analysis, as the number of degrees of freedom increases, can 
be found in Figure [H] and Table O Due to the discussion of Section 14. 4^ 
we consider this for different ratios N/M. 

The right graph shows that the conclusion of Corollary [6] is far from 
optimal for the inverse harmonic potential of Section 14. 3[ As expected 
from Section 14.41 a faster convergence rate as well as smaller residuals 
are found if we suitably choose < M. 

The left graph corresponds to the coulombic potential discussed in 
Section 14.11 It clearly shows that the order of convergence of A„ — * 
does not obey the estimate |A„ — E\ < 0{n~^) (for some a > 0) of 
Corollary |H In fact the convergence rate seems to decrease as we in- 
crease the number of degrees of freedom. This reduction in the speed of 
convergence can be prevented by putting M = f{N) for a suitable non- 
linear increasing function < f{x) < x. The optimal f{x), however, 
might depend on the eigenvalue to be approximated. 

Remark 6. According to Remark [2], the actual approximate eigenvalue 
Re(A) is correct up to 0(?7,~^"), where a is the second column of TableO 
Furthermore, note that, in the case of the coulombic potential we can 
compute directly the true residual | Re(A) — E\. From Figure [HI bottom, 
it is clear that this true residual is substantially smaller than the one 
estimated by | Im(A)|. 
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Appendix A. Entries of the matrix polynomial 

COEFFICIENTS 

The recursive identities satisfied by the Hermite functions allow us 
to find recursive expressions for the matrix entries of K and L in (13.11) 
when G = H^. Rather than estimating the corresponding inner prod- 
ucts by trapezoidal rules, we build the codes involved in the numerical 
experiments performed in Section H] using these explicit expressions. As 
large factors are cancelled in these explicit expressions, this approach 
turns out to be far more accurate. Since some of the calculations are 
not entirely trivial, we include here the crucial details. 

Let 

/•oo poo 

Ti= / <l>,.(r)<l>,(r)dr, T2[k,])= $'fc(r)$,(r) dr, 
JO Jo 



Ts 

Fi -- 



-(^i.{r)^Ar) dr. 



r 

oo - 



-<l>fc(r)<l>j(r) dr. 



V)$fc(r)$j(r) dr. 



T4 = 
F2 = 



$fc(r)$; (r) dr. 



oo -j^ 

r 



— $fc(r)<l>j (r) dr. 



6^i(r)$fc(r)$j(r)dr, 



$fc(r)$j(r)dr, Fi{k,j) 



5el 



V)$'fc(r)$j(r) dr. 



Here and below we stress the dependence on j, k when the coefficient 
is not symmetric with respect to these indices. Denote 







j,2 



Then {H^-^u, ^jm) are given according to Table [T] and {H^'^ki, -f^K^jm) 
are given according to Table [2l 





m = 1 


m = 2 


/ = 1 


Ti + F, 


T2{k,j) + kTs 


/ = 2 


-T2{k,j) + kTs 


-Ti + Fi 



Table 1. Term {H^^kh^jm)- 



For m, n e N U 0, let 



Pin) 



nr=i(i + ^) n^O 

1 n = 
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m = 1 


m = 2 


/ = 1 


Ti + T4 + Kmk,j) + kt^u, k)+ 

K^Tg + 2Fi + F2 


-T2ik,j)-T2ij,k) + 2KFs+ 
F,ik,j)-F,{j, k) 


/ = 2 


-T2ik,j)-T2ij, k) + 2KFs+ 
~F,ik,j) + F^{j,k) 


T1 + T4- nT^^kJ) - nT^ij, k) + 
K^Te - 2Fi + F2 



Table 2. Term {H^^kuH^^jrn)- 



and 



(A.l) I{m,n) = / hm{r)K{r)e dr. 



Cm C' 



m^n Jo 



Lemma 7. 



6mn m = n (mod 2) 



(2fc-2j-l)^7r(2fc+l) 

Proof. \im = n (mod 2), then hm{r)hn{r) is an even function for r G 
and so 



1 







hm{r)hnir)e dr = - / h^{r)hn{r)e dr. 



00 



On the other hand, if m ^ n (mod 2), say m = 2k and = 2j + 1, 
f l3.10p and integration by parts yield 



h2k{r)h2,+i{r)e-'" dr = / /i2,+i(r)(e-'-')(2'=) dr 
Vo 

00 

2 N 







/■oo 

-2(2j + l) / /i2j(r)(e-'-')(2'=-i)dr 

/>oo 

22(j + l){2k - 1) / /i2fc-2(r)/i2,-i(r)e-'^' dr. 



The corresponding expression for I{m, n) can be obtained in a straight- 
forward manner from these two assertions. □ 
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Lemma 8. 

Ti(k,j) = 6,k, T2(k,j) = ^, '^y--^') -JpifsPik), 

^ ' ^ 0r(2A;-2j-l)(2A:-2j + l)^ ^ ^' 



{-v/2A;(2fc + 1) J = A;-1 

4^ + 3 J = fc 

-v/(2A; + 2)(2A; + 3) j = A; + 1 

otherwise. 



r5(fc,j) = <; 1 k = 3 T,{k,j) = i-iy-'2 

A; > j, I V ^ 

Proof. Let /(m, n) be given by (lA.ip . By virtue of identities (13.101) and 

POO 

/ Kk+ih2j+ie-'' dr = {2k + l)I{2k, 2j + 1) - 2-'l{2k + 2, 2j + 1), 
Jo 

- W/i2,+ie-^' dr = ^(-1)'22'+^— --/(2(A: - /), 2j + 1). 

1=0 \ )■ 

This yields T2 and T3. 
Let 

J{k,j) = / -h2kh2j+ie dr. 
Jo 

Then 

f 0F2^J(2fc)!(-l)J~fej! 



k> J 



and 

r"00 



POO 

/ -h2k+ih2j+ie-'' dr = 2(2A: + 1) J(A;, j) - 4iv^2*^(2fc + 1)!. 
Jo 

This renders T5. Moreover, integration by parts ensures 

n = I = - (-) $ = {T,{k,j) + uj, k)). 

The expression for T4 follows from fl3.13l) and the identity 



$'fc(r)$;.(r) dr = - / -$;'(r)<l>j(r) dr. 

□ 
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From El and E2 in the next lemma, one easily obtains explicit for- 
mulae for Fn when </'ci(r) = '-f/r^. 

Lemma 9. For (3 G [0, 1] and a G [-1, 2], let 

POO POO 2 

Eli/3, k,j) = / ^<l>l(r)%(r)dr, E2{a,k,j)= / -$fc(r)<l>,(r) dr. 
Jo Jo ^ 

Then 



Ei{P,k,j) = {2j+l)E2{P+l,k,j) + ^/2{2j + l)jE2{P+l,k,j-l)-E2{P-l,k,j) 
and 



^ m\p\P(m)P{p) 

Proof. Let 



= C2,\A2k + 1)1- 



' {k - n)\{2n + 
Then 



C2,'+M+i{r) = J2Sn{ky-^' 



n=0 

and 



^2(«,fc,j)= ^ Smik)Sp{j)K{a,m,p) 

m,p=l 



where 



poo I 1 / Q n \ 

K{a, m,p)= / —r^rn+i^^p+i^-r^ ^ -r ( — — + m + p) . 

Jq 2 V 2 / 

On the other hand, the expression for Ei follows from applying (13.101) 
and fl3TT]) . □ 
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If and E4 are as in the following lemma and 0ci(^) = then 
F,{k,j) = Es{2k + l,2j + l), 

1 / l2k + 1 l2i + 1 

F2{K3) = - (^^ ^—E,i2k,2j + 1) + ^-L-—E,i2k + l,2j) 



- VkTTEi{2k + 2, 2j + 1) - v^J + 1^4 (2 A; + 1, 2j + 2)- 
E3(2fc + l,2j + l)' 



/ 2/;; + 1 

F3(A:,j) = T,ik,j) - J^—E3i2k,2j + 1) - v/^E3(2A; + 2, 2j + 1), 



/2A; + 1 

F^ik, 3) = y -^Es{2k, 2j + 1) - v/^E3(2A: + 2, 2j + 1). 
Lemma 10. For m^n G NU {0}, /ei I{m,n) be as in Lemma^ 

1 r 1 



E:i(m,n) 



E^{m, n) 



CmCn Jo 1 + T 



CmCn Jo ^ ^ 



;hmir)hn{r)e dr, 
jhm{r)hn{r)e~''^ dr. 



Then 



E3(0,0) = 2e/" e-*'dt + ev^, 
Jo 

^4(0,0) =ev^y — dt, 



^3(m + l,0) 
E4(m + 1,0) 
^3(771 + l,n + 1) = 

£■4(771 + 1,77 + 1) = 



^i== (£4(771, 0) - V^Esim - 1, 0)), 
y/m + 1 

^ ^ f 72/(771, 0) - V2E3{m, 0) - y/mE^im - 1, 0)) 

V777 + 1 V / 
: (2/(777, 77) - 2£'3(777,77) + 



v/(7r7+ 1)(77+ 1) 

— V 2mE4{m — l,n) — \/2nE4^{m, 77 — 1) + y/mnE^lm — l,n — 1)) 
— =^=== (v/2(72 + l)/(m, 77 + 1) - 2£4(m, 77) + 

-y/(7?7 + 1)(77 + 1) ^ 

^2?7/?3(m, 77 — 1) — V 2777/^4(777 — 1, 77) + \/ ninE3{m — 1, — 1) j . 



Proof. The recursions for £^3 and /?4, follow from (13.111) . 



□ 
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Figure 1. Portion of Spec2(i/_i, in thin boxes 
around [—3, 3] for the purely coulombic potential with 
7 = —1/2. The bottom image shows details of the 
picture near 1. Here we superimpose two values of 
n, 500 and 1000. According to Theorem [T] there are 
approximate energy states at ~ 0.8661 ± 0.0236, 
E ^ 0.9662 ± 0.0086 and E ^ 0.9853 ± 0.0041. These 
correspond to the actual eigenvalues Eo ^ 0.866025, 
El ^ 0.965925 and E2 ^ 0.985121. 
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Figure 2. Approximate ground wave function for a 
purely coulombic potential with 7 = —1/2. The true 
wave function (blue line) can be found explicitly, |Tha92l 
Section 7.4.2]. The ground state in this case is Eq = 
^1 — 7^ ^ 0.866025. We have deliberately chosen small 
dimensions n for the test spaces, in order to illustrate ap- 
proximation of the method. The residual error is actually 
much smaller that the one predicted by Theorem [H see 
Table O 

^ \\v-Usv\\ |Im(A)| 
dE 

15 0.176115 0.680599 

25 0.084527 0.514205 

35 0.072552 0.457034 
Table 3. Here we compare both sides of (13.41) . for the 
computation of the approximate eigenfunctions of Fig- 
ure El We approximate (Ie = Ei - Eq ^ 0.0999004. 
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Figure 3. Computation of the ground energy value for 
and (l)e\{r) = We depict Eq against p. The 

vertical bars correspond to the error predicted by Theo- 
rem [H 
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Table 4. Data depicted in Figure [3] 
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Figure 4. The graph captures the evolution of as it 
crosses the spectral gap of where (\)c\{t) — 7/(1 + 
r^) for 7 = —5 : .5 : 0. We consider three choices of 
pairs {N^M) such that d\m{CMM) = 120. The curve 
corresponds to Re(An) for A = A„ in (13 ■4p and (13.151) . 
The vertical bars on the curve measure |Im(A„)| . We 
superimpose the image with the eigenvalues of L in (13.11) 
for G = H_i, that is the Galerkin approximation. 
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Figure 5. Here we use Theorem [T] to find the first three 
eigenfunctions of with 4>ei{f) = —4/ (1 + r^). The nu- 
merical evidence suggests: Eq ^ —0.3955, Ei ^ 0.6049 
and E2 f» 0.9328. 
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Figure 6. Here Eq and Ei are eigenvalues of for 
(pei{r) = — ^. The top graph shows the eigenvalues of 
L in (13. 1!) (that is the Galerkin approximation) for G = 
H^i and (M, iV) = (iV, 200 - N) so that dim{CNM) = 
200. The bottom graph depicts the residuals |Im(A„)| 
and I Re(A) — Ej\. For Eq, the minimum of the residual 
curve corresponding to |Im(A)| is achieved when ^ 
155 and it is roughly 7% smaller than when = 100. For 
the same eigenvalue, the residual curve corresponding to 
I Re(A) — Eo\ achieves its minimum when = 165 and 
it is roughly 66% smaller than when N = 100. 
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Figure 7. Here E'o, Ei and E2, are the first three eigen- 
values of H^i for 0ei(^) = ~4/(l + r^). The top graph 
shows approximation of Eq ^ —0.3955, Ei ^ 0.6049 
and E2 ^ 0.9328, for (M, N) = {N, 120 - A^) so that 
dim(£7VM) = 120. The curves correspond to Re(A„) for 
A = A„ in (13. 4p and (13.151) . The vertical bars measure 
|Im(A„)|. The image is superimposed with the eigen- 
values of L in (13. ip for G = H_i, that is the Galerkin 
approximation. The bottom graph depicts the residuals 
|Im(A„)|. 
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Figure 8. Log-log plots of |Im(A)| for Re(A) close to 
an eigenvalue, Eq, for different choices of pairs (A^, M) 
a.s n = N + M increases. Left: k = —1, (pcii'f') = 
and Eo ^ 0.86602. Right: k = -1, </)ei(r) = --^ and 
Eq ^ 0.61399. See Table El 
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Table 


5. In this table we 


fit by least squares the dat 



of Figure M and find a and b such that |A„ — £'o| < 
I Im(A„)| ~ bn" ioT n = N + M. 
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